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1  Summary 


The  goals  of  this  project  ha%e  been  defined  in  1985  as  follows:  “An  improved  global  cali¬ 
bration  will  be  derived  and  image  enhancement  techniques  developed  to  increase  the  image 
quality,  sensitivity,  and  spatial  resolution  of  maps  from  the  Infrared  Astronomy  Satellite 
(IRAS)  survey  data.  Enhanced  maps  of  selected  regions  of  the  sky  will  be  used  to  improve 
and  validate  the  AFGL  infrared  celestial  background  model”.  These  goals  were  achieved  by 
augmenting  and  extending  the  work  started  in  1984  on  IRAS  data  Processing  in  Groningen. 
The  total  project,  operational  from  1984  to  probably  1992,  is  called  “GEISHA”  (Gronin¬ 
gen  Exportable  IRAS  High-resolution  Analysis)  system.  All  goals  of  the  AFOSR  part  of 
GEISHA  were  achieved  by  March  '89.  Therefore,  it  was  necessary  to  make  a  few  more 
products  than  anticipated. 

A  completely  new  calibration  scheme  using  the  zodiacal  light  emission  as  a  standard 
candle  has  been  designed  and  implemented.  The  standard  IRAS  data  products  presently 
available,  that  were  produced  at  the  IR.A.S  Processing  and  Analysis  Center  (IPAC  at  Caltech. 
Pasadena,  U  S. A  ),  were  all  calibrated  using  the  point-source-like  response  of  the  calibration 
“flashes”.  Our  calibration  is  better  suited  for  smaller  spatial  frequencies.  Comparison  of  the 
zodiacal  light  calibration  and  that  of  IP.-\C.  for  point  sources,  gives  very  satisfactory  results 
(<  10  %)  at  the  wavelengths  12,  25,  and  60  fim,  but  a  50  %  error  at  100  ^m.  Also  scans 
making  an  angle  larger  than  5“  with  the  plane  perpendicular  to  the  line  Sun-IRAS  cannot 
be  calibrated  properly;  there  are  only  few  of  such  scans. 

An  iterative  least  squares  (ILSQ)  routine  has  been  prepared  that  makes  it  possible  to 
enhance  the  resolution  of  an  IRAS  image  to  slightly  better  than  2  arc  minutes,  for  most 
regions  of  sky.  The  output  of  this  routine  is  directly  traceable  to  the  input  and  the  method 
is  rather  straightforward  and,  because  of  that,  a  lot  faster  than  other  non-linear  methods 
developed  to  restore  IRAS  images;  e.g.  FIER  by  Kennealy  et  al.  (1987).  FIER  produces 
higher-resolution  images  than  ILSQ,  but  the  method  needs  verification.  Our  ILSQ  program 
is  ideal  for  creating  large  spatial  resolution  enhanced  images,  and  for  checking  e.g.  FIER 
results. 

While  preparing  these  two  main  products  it  became  evident  that  a  number  of  other  — 
originally  unforeseen  —  products  were  also  needed. 

It  was  not  as  simple  as  expected  to  combine  the  detector  data  and  the  pointing  restoration. 
The  IPAC-produced  Boresight  Pointing  History  File  (BPHF)  needed  to  be  reduced  by 
a  factor  five,  and  to  be  cut  in  pieces  that  corresponded  to  the  363  sky  regions  defined 
in  an  earlier  stage  of  the  GEISHA  project  (stored  on  363  “plate”  tapes). 

The  “sample”  file  was  defined,  a  combination  in  one  Fits-like  file  of  the  calibrated  detector 
data  and  the  pointing  information. 

Before  being  able  to  run  the  ILSQ  program  it  is  essential  to  correct  the  input  image  for 
responses  not  caused  by  actual  sky  sources  (deglitching)  and  for  the  baseline  differences 
between  scans  (destriping).  A  deglitching/destriping/imaging  program  was  made,  and 
is  so  succesful  that  the  IPAC-produced  Super  Sky  Flux  product  uses  that  routine. 

A  spline-smewthed  “plate”  system  was  created  as  well,  consisting  of  40  instead  of  363 
tapes.  Images  with  a  resolution  of  8  arc  minutes  can  be  made  using  this  data  base, 
and  because  of  the  low  spatial  resolution  it  is  feasible  to  make  images  of  very  large  sky- 
regions. 
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Resolution-enhanced  images  of  quite  a  number  of  sky  regions  have  been  produced  already, 
and  are  being  produced  at  present.  The  recent  arrival  of  an  Alliant  computer  will  make  it 
possible  to  handle  really  large  regions  of  sky. 

An  overview  of  the  GEISHA  system  has  been  published  in  1988  (see  Appendix  1).  It 
now  runs  on  a  Tybpr  170/760  (NOS),  a  VAX8600  (VMS),  a  SUN  3-280  (UNIX),  an  Alliant 
(UNIX),  and  a  Convex  (UNIX).  This  proves  lhai  ihe  package  is  exportable,  although  it  is 
advisable  —  at  present  —  to  use  sample  files,  created  in  Groningen,  as  input  for  further 
processing. 

2  Data  Bases 

The  following  data  bases  have  been  created  and  will  be  maintained  in  the  coming  years. 

2.1  complete  scans  in  FITS  format 

We  originally  started  the  job  with  1000  magnetic  tapes  imported  to  Groningen  from  the 
IRAS  ground  station  at  Chilton,  U.K.  in  1984.  As  the  first  data  processing  step  an  index  of 
the  contents  of  these  tapes  was  made  and  they  were  converted  to  the  astronomical  standard 
for  data  exchange.  FITS  (Flexible  Image  Transport  Standard).  That  led  to  6  sets  of  30 
tapes  for  the  up  and  down  scans  of  the  three  HCON’s,  2  sets  of  5  tapes  for  the  Minisurvey. 
60  tapes  for  the  AO  data,  and  a  few  tapes  for  calibration  and  administration  (these  have 
the  labels  RID<id><nn>;  id  =  A,  B,  C,  D,  E,  F  for  the  survey  scans;  id  =  G.  H  for  the 
minisurvey:  id  =  I  for  AO's;  id  =  J  for  calibration  and  administration). 

2.2  Plate  Tapes 

These  complete  scans  were  subsequently  cut  into  pieces  of  about  15°  long  and  sorted  to  sky 
regions.  The  goal  of  the  sorting  operation  was  to  put  all  IRAS  data  pertaining  to  a  certain 
region  of  sky  (of  order  15°  X  15°)  on  one  magnetic  tape.  Another  requirement  w  as  to  be 
able  to  make  a  map  of  any  1°  x  1°  region  of  sky  just  accessing  one  tape.  In  order  to  achieve 
this  we  had  to  divide  the  sky  in  363  regions,  and  consequently  the  data  base  consists  of  3C3 
magnetic  tapes  (density  6250  bpi). 

If  a  sky  region  needs  to  be  studied  that  is  larger  than  1°  x  1°  it  may  be  needed  to  access 
more  than  one  tape.  For  convenience  in  transporting  the  data  of  such  regions  there  exist 
cutting  and  glueing  programs  to  create  a  new  tape  (or  set  of  tapes)  that  contain  all  data  of 
the  desired  sky  region. 

The  actual  sorting  to  plate  tapes  was  done  in  two  stages,  because  it  was  more  efficient 
and  to  have  another  intermediate  data  base  as  a  backup  (tapes  are  definitely  not  the  best 
way  to  store  data  permanently).  In  the  first  stage  the  scans  were  cut  into  scan  snips  suitable 
to  incorporate  in  the  plate  tapes,  but  only  ordered  to  families,  i.e.  rather  large  sky  regions. 
That  led  to  tapes  RIF<Fam.><nn>  (Fam.  =  A  . .  .S). 

Subsequently  each  large  sky  region  was  sorted  to  plates,  the  tapes  RIPOOl  . .  .RIP363. 

2.3  BPHF 

The  boresight  pointing  history  file  (BPHF),  giving  information  about  the  pointing  direction 
of  IRAS  for  every  second  of  its  life  (about  26  million  seconds)  has  been  obtained  from  IPAC, 
103  magnetic  tapes  (BPHOl  ...BPH103).  By  omitting  redundant  numbers  and  recording 
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movement  rather  than  absolute  position,  it  was  possible  to  reduce  this  data  base  to  21  tapes 
in  FITS  format  (RIBHOl  . .  .RIBH21,  at  1600  bpi),  without  loss  of  accuracy. 

The  BPHF  was  cut  into  snips  and  sorted  such  that  all  BPHF  data  relevant  to  a  certain 
plate  were  combined.  It  is  the  eventual  aim  to  put  these  data  on  the  plate  tapes  themselves, 
but  at  present  they  are  stored  on  separate  tapes  (RIP401  .  ..RIP408). 

Many  sky  maps  have  already  been  made  by  combining  detector  data  and  BPHF  data,  and 
usually  that  is  quite  succesful.  We  have  verified  that  the  match  is  normally  OK.  However, 
about  one  scan  per  region  of  3°  X  3“  lies  of  the  order  of  15°  away  from  the  region.  We  have 
been  unable  until  now  to  discover  where  the  error  is  introduced.  In  the  near  future  we  will 
obtain  SUPERBORESIGHT  from  IPAC.  Then  we  will  snip  and  order  those  data  anew,  and 
hopefully  the  problem  will  disappear. 

2.4  Miscellaneous 

Some  other  data  bases,  relevant  for  IRAS  data  processing,  are  stored  bls  well. 

The  response  functions  of  all  detectors  as  derived  by  Moshir  from  IPAC  are  on  tape 
RRBOOl. 

All  raw  detector  data  were  smoothed  using  cubic  splines  to  a  resolution  of  about  8  arc 
minutes.  Then  all  IRAS  data  can  be  put  on  just  25  tapes,  RIDS<Hcon><n>. 

Finally,  the  data  pertaining  to  all  stimulator  flashes,  that  were  done  at  the  begin  and 
the  end  of  each  scan,  were  put  aside  and  compressed,  by  using  an  appropriate  template,  to 
just  one  value  and  its  uncertainly  per  flash.  This  led  to  two  tapes  (RIFZCO,  RIFZCl). 

3  Sample  files 

In  order  to  make  the  raw  IRAS  data  in  a  computer-independent  way  available  to  the  astro¬ 
nomical  community  a  new  intermediate  data  product  has  been  designed:  sample  files.  For 
a  given  region  of  the  sky,  a  sample  file  contains  all  information  present  for  one  of  the  IRAS 
wavelength  bands.  Each  flux  calibrated  detector  readout  and  its  position  are  present, 
leading  to  a  long  table  of  brightnesses  and  coordinates,  and  their  respective  uncertainties. 
Even  the  orientation  on  the  sky  of  a  particular  detector  is  given  for  each  sample. 

This  table  can  be  processed  further  for  several  applications.  The  main  one  is  a  destrip- 
ing/imaging  program,  followed  by  the  iterative  least  squares  image  enhancement  program. 

It  is  also  possible  to  study  individual  scans  by  plotting  intensity  versus  time,  and  to 
coadd  the  data  of  all  scans  running  through  a  particular  point  on  the  sky. 

4  Destriping  and  Imaging 

The  IRAS  satellite  made  scans  of  order  180°  long  and  0.5°  wide.  The  first  set  of  maps  were 
made  at  IPAC,  the  Sky  Fluz  maps.  These  have  been  widely  used.  They  were  made  by 
dividing  the  sky  in  bins  of  2  x  2  arc  minutes;  when  a  sample  overlapped  a  bin  its  value 
was  stored  in  the  bin.  Averaging  all  values  in  a  bin  resulted  in  a  sky  map.  Maps  created  in 
this  way  have  a  resolution  corresponding  to  the  detector  size  or  poorer,  if  scams  at  different 
angles  are  combined.  Other  methods  to  co-add  the  scans  lead  to  similar  resolutions. 

Imaging  methods  like  those  creating  Sky  Flux  leave  responses  due  to  stripes,  spikes, 
and  space  debris  in  the  map.  Stripes  are  caused  both  by  local  imperfections  in  the  flux 
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calibration  which  show  as  straight-lined  structure  in  the  scan  direction  of  the  satellite  and 
by  differences  in  the  zodiacal  light  level  at  the  time  of  observation. 

Spikes  are  the  results  of  particle  hits  on  the  detectors.  They  can  be  mistaken  for  genuine 
point  sources  (a  point  source  is  only  seen  in  two  to  three  samples).  Not  too  close-by  space 
debris  and  asteroids  are  indistinguishable  from  point  sources  in  a  single  scan.  The  IRAS 
Point  Source  Catalog  had  a  confirmation  strategy  that  eliminated  all  such  undesired  point 
sources,  but  for  images  other  methods  have  to  be  used. 

A  program  has  been  developed  to  produce  images  from  sample  files  and  to  correct  these 
images  in  an  interactive  way  for  striping  and  unconfirmed  point  sources.  It  is  possible  to 
read  in  an  existing  map  (in  FITS  format)  or  the  program  can  create  a  new  map  first.  Either 
just  a  map  or  an  iteratively  destriped  map  can  be  made.  Images  are  produced  by  a  coadd 
of  sample  file  data  into  a  map  either  with  a  circular  beam  or  with  the  ’true’  detector  size 
as  a  beam.  Some  scans  can  be  excluded  from  the  process,  if  they  are  deemed  bad  or  are 
unwanted  otherwise. 

Stripes  are  removed  as  follows.  Each  individual  detector  scan  is  modelled  as  a  straight 
line.  The  line,  determined  using  linear  regression,  is  subtracted  from  the  detector  scan  and 
the  result  is  coadded  into  a  new  map.  This  process  is  done  for  all  scans.  A  new  map  ensues 
which  contains  less  stripes.  Typically,  the  whole  procedure  is  iterated  five  times.  Usually 
the  resulting  maps  have  no  visible  stripes  any  more.  Occasionally  the  map  is  still  bad  in  a 
small  region.  That  is  often  due  to  calibration  problems  with  the  snip  of  an  originally  short 
scan  (not  running  from  ecliptic  pole  to  pole).  A  simple  remedy  is  removing  that  scan  snip. 

Spikes,  asteroids  and  space  debris  are  eliminated  from  the  final  map  by  applying  a 
median  filter  in  the  coadd  algorithm.  The  unwanted  responses  appear  only  in  some  of  the 
overlapping  scans,  while  the  majority  of  the  scans  has  proper  information  on  that  position. 
A  median  filter  puts  heavy  emphasis  on  the  response  of  the  majority  of  the  scans,  and 
eliminates  largely  deviating  points. 

With  this  program  it  is  possible  to  make  sky  maps  retaining  the  resolution  of  a  simple 
coadd.  The  maps  are  free  from  artefacts  which  plagued  the  Sky  Flux  maps  such  as  stripes, 
spikes,  space  debris,  and  asteroids.  Moreover,  all  three  HCONs  can  be  combined  into  one 
map  which  increases  the  signal-to-noise  ratio  by  at  least  a  factor  2.  Finally,  the  output  from 
this  program  is  the  ideal  input  to  the  ILSQ  program  described  below. 

5  Calibrating  IRAS  using  Zodiacal  Emission 

Inspired  by  the  work  of  Dr.  Price  during  a  visit  of  half  a  year  to  Groningen  a  new  calibration 
has  been  developed  employing  the  zodiacal  emission,  that  is  clearly  present  in  each  long  scan 
(at  least  at  12,  25  ,  and  60  /rm).  Such  a  calibration  has  advantages  over  the  scheme  employed 
by  IPAC.  The  IPAC  calibration  uses  the  very  short  (few  tenth  of  a  second)  stimulator  flashes 
at  the  begin  and  end  of  each  scan.  These  are  ideal  to  calibrate  the  detector  response  to  point 
sources,  but  it  is  known  that  the  long-term  behaviour  of  the  detectors  is  different.  In  order 
to  calibrate  extended  features  using  the  zodiacal  emission  that  extends  over  many  tenths 
of  degrees  as  a  stamdard  candle  might  produce  a  more  consistent  result.  IPAC  applied  the 
point  source  calibration  to  their  Sky  Flux  product.  It  is  interesting  to  compare  our  eventual 
maps  with  theirs. 

We  have  assumed,  because  of  several  IRAS  studies  of  zodiacal  emission  by  others,  that  the 
zodiacal  emission  as  a  whole  can  be  approximated  by  the  emission  of  a  tilted  oblate  ellipsoid 
of  uniform  density  and  temperature.  The  socalled  side  bands  at  low  ecliptic  latitudes  have 
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newly  been  discovered  by  IRAS,  are  interesting  in  their  own  right,  but  are  a  nuisance  for 
our  calibration  goal.  We  therefore  ignore  a  small  part  of  the  emission  near  ecliptic  latitude 
0“. 

Apart  from  the  zodiacal  emission,  there  is  always  a  contribution  of  gadactic  emission, 
especially  prominent  at  60  and  100  pm.  The  large-scale  galactic  emission  can,  just  as  the 
zodiacal  emission,  be  used  as  a  standard  candle.  At  100  pm  the  large-scale  galactic  emission 
is  in  fact  more  important  for  calibration  purposes  than  the  zodiacal  emission,  which  is  faint 
and  only  present  at  low  ecliptic  latitudes  (just  a  few  degrees).  The  method  is  described  in 
more  detail  in  a  conference  paper  that  is  attached  as  Appendix  2. 

To  calculate  the  calibration  parameters  we  make  the  following  assumptions.  As  a  stan¬ 
dard  candle  we  take  the  surri  of  a  (too  be  determined)  Zodiacal  Emission  Model  (ZEM)  and 
Galactic  Emission  Model  (GEM):  ZEM  -t-  GEM .  The  baseline  is  allowed  to  change  linearly 
over  a  scan.  It  is  assumed  that  the  slow  exponential  decay  of  the  detector  response  that 
occurs  will  be  sufficiently  approximated  by  a  straight  line.  Finally,  an  in-house  developed 
model  for  the  non-linear  behaviour  of  the  gain  is  applied. 

All  linearization  corrections  to  the  data  have  to  be  applied  first.  The  actual  model  values 
for  ZEM  and  GEM  have  to  be  derived  from  the  linearized  IRAS  uncalibrated  data  in  a  kind 
of  least-squares  fashion.  Many  iterations  through  all  the  data  of  a  band  are  needed  to  derive 
a  final  calibration. 

To  actually  perform  this  calibration  a  lot  of  data  processing  had  to  be  done,  and  software 
had  to  be  created.  The  original  uncalibrated  IRAS  data  were  smoothed,  in  a  number  of  steps, 
to  just  a  few  promiUc  of  the  data  volume.  These  smoothed  data  were  extensively  studied  to 
arrive  at  a  suitable  ZEM  and  GEM,  and  to  develop  a  processing  method. 

The  smoothed  data  base  still  comprises  several  hundred  Mbytes,  and  running  one  cal¬ 
ibration  task  for  one  band  (a  lot  of  iterative  operations  are  done)  costs  a  few  days  on  a 
stand-alone  SUN  3-280.  Therefore  a  simplified  version,  only  using  a  ZEM  as  the  standard 
candle,  is  operational  at  present. 

The  simple  calibration  compares  well  with  the  Point  Source  Catsdog  calibration  at  12. 
25,  and  60  fxm  (to  within  10  %).  It  does  not  work  well  if  the  scans  to  be  calibrated  had  an 
angle  to  the  Sun-satellite  direction  that  deviates  more  than  ten  degrees  from  90°.  It  also 
does  not  work  well  at  100  /ira,  because  the  zodiacal  emission  is  a  feeble  standard  candle  at 
that  wavelength. 

In  the  summer  months  of  1989  a  much  larger,  AUiant  computer  will  become  available 
and  will  make  it  possible  to  run  the  calibration  jobs  in  a  reasonable  time.  We  expect  to  have 
the  improved  ZEM  -1-  GEM  calibration  ready  and  operational  by  October  1989. 

6  ILSQ  Image  Enhancement  Program 

Reconstruction  of  images  from  a  collection  of  scans,  like  obtained  with  IRAS,  requires  the  de- 
velopement  of  new  techniques  in  the  field  of  image  restoration.  Conventional  image  restora¬ 
tion  methods  do  not  work,  because  they  start  from  already  existing  images,  which  are 
degraded  to  some  extent.  A  collection  of  scans,  however,  does  not  form  an  image.  A  scan 
can  be  regarded  as  a  one-dimensional  cross-cut  through  the  image.  More  precisely,  a  scan 
represents  the  brightness  of  a  narrow  strip  of  the  sky,  convolved  with  the  response  function 
of  one  of  the  IRAS  detectors.  The  image  has  to  be  reconstructed  from  the  deconvolution  of 
the  scans,  together  with  the  arrangement  of  the  scans  over  the  region. 
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It  is  especiall>  important  to  improve  the  poor  spatial  resolution  in  the  cross-scan  direc¬ 
tion.  That  is  —  in  principle  —  possible,  because  of  the  confirmation  IRAS  survey  strategy; 
over  most  positions  on  the  sky  there  run  six  scans,  and  more  specifically  the  detectors  are 
shifted  by  half  a  detector  width  between  two  —  one  orbit  apart  —  scans. 

An  iterative  bnear  least  squares  deconvolution  scheme  has  been  devised  which  is  insensi¬ 
tive  to  the  uneven  coverage  of  an  area  by  scans.  The  method  allows  both  to  obtain  a  higher 
angular  resolution  than  the  simple  co-adding  technique  and  to  provide  an  ensuing  map  on 
a  regular  mesh.  In  order  to  apply  this  method  without  introducing  artefacts  it  is  needed  to 
use  the  best  possible  response  functions  of  the  detectors;  the  rectangular  approximation  of 
the  co-adding  ("Sky  Flux”  type)  method  is  inadequate.  The  method  makes  the  following 
assumptions. 

Each  sample  is  regarded  to  be  the  convolution  of  the  sky  with  the  response  function  of 
the  detector  that  has  obtained  the  sample.  It  is  assumed  that  the  zodiacal  emission  has 
already  been  subtracted,  and  that  the  detection  process  is  linear,  which  ignores  the  memory 
effects  in  the  detectors.  The  map  area  is  divided  into  a  number  of  pixels,  and  the  brightness 
is  regarded  constant  over  one  pixel  area.  The  pixels  are  surfaces  and  measure  typically  onf" 
square  arc  minute. 

The  ILSQ  method  amounts  to  solving  an  equation  of  the  type: 

D  =  PB,  (1) 

with  D  a  vector  of  length  A’  representing  all  measured  values,  B  a  vector  of  length  M  being 
the  unknown  brightnesses,  and  P  i  N  x  A/  matrix.  Clearly  for  this  equation  to  be  solvable 
N  >  M;  also  all  pixels  must  be  covered  at  least  once. 

The  possible  angular  resolution  of  the  reconstructed  image  depends  on  the  number  of 
samples  per  unit  area.  The  plac"  and  orientation  of  the  scans  running  over  a  map  with 
respect  to  each  other  will  obviously  also  play  a  role;  the  sought  spatial  frequencies  must  be 
present. 

Fortunately  the  matrix  P  is  sparse,  because  the  solid  angle  of  the  response  function  of  one 
detector  covers  an  area  much  smaller  than  the  usuadiy  requested  map  area.  The  majority 
of  pixels  is  not  covered  when  a  sample  is  overlayed  on  the  map  grid.  The  sparseness  is 
determined  by  the  surface  area  ratio. 

Direct  least  squares  methods,  which  solve  overdetermined  systems,  tend  to  fill  the  entire 
matrix  P,  which  becomes  rapidly  unmanageable  on  present  computers.  Iterative  methods 
that  preserve  the  sparseness  of  the  matrix  are  therefore  required.  The  LSQR  algorithm  of 
Paige  and  Saunders  (1982)  turns  out  to  work  well,  and  does  even  more,  viz.  it  solves  damped 
least  squares  problems.  The  introduction  of  a  damping  factor  A  ‘regularises'  the  problem, 
in  the  sense  that  the  solution  B  has  become  less  sensitive  to  outliers  in  the  data  D. 

The  right  damping  factor  hais  to  be  chosen  for  each  map  anew.  The  pixel  size  should  be 
chosen  such  that  the  system  remains  sufficiently  overdetermined:  the  smaller  the  pixels  the 
larger  the  noise.  In  general,  the  reconstruction  of  a  map  turns  out  to  be  acceptable  when 
the  number  of  equations  N  exceeds  the  number  of  pixels  A/  by  a  factor  3-4. 

We  conclude  that  the  iterative  LSQ  method  to  produce  images  of  the  far-infrared  sky 
using  the  IRAS  survey  data  as  input  works  very  satisfactorily.  For  most  regions  of  sky, 
images  can  be  made  with  a  resolution  of  slightly  better  than  2  arc  minutes,  also  in  the 
cross-scan  direction.  A  more  extensive  paper  on  this  topic  is  in  preparation  and  some 
preliminary  figures  of  that  paper  illustrating  the  results  have  been  attached  to  this  report 
as  Appendix  3. 
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SUMMARY 


The  GEISHA  (Groningen  Exportable  Infrared  Survey  Hlgh-resolution  Analysis)  sys¬ 
tem  has  as  its  two  main  goals  to  provide  astronomers  with  direct  access  to  the 
raw  (uninterpreted)  data  and  to  develop  high-spatial  resolution  algorithms.  Ex- 
portability  was  formulated  as  a  secondary  design  goal  of  the  prcjcc*’. 

Tne  problem  of  calibration  of  detector  read-outs  in  terms  of  surface  brightness 
distribution  over  tne  sky  involves  a  non-trivial,  and  possibly  subjective, 
decomposition  in  detector  physics  and  astrophysics.  This  led  us  to  the  decision 
to  make  un-lnterpreted  detector  readouts  accessible,  especially  for  attempting 
the  highest  possible  resolution  in  Imaging,  but  also  for  study  of  large-scale 
structures  (tens  of  degrees). 

GEISHA  has  two  main  output  products  for  astronomers:  sample  files  and  sky  maps. 
A  sample  file  has  ail  necessary  information  on  detector  positions  and  orienta¬ 
tions,  errors,  and  flux  calibration.  It  can  be  used  to  sake  a  sky  map,  or  as 
input  for  the  application  of  high-'resolution  techniques. 

In  order  to  be  able  to  make  maps,  surface  brightness  and  position  calibrations 
are  required.  If  (extra)galactlc  emission  needs  to  be  studied  a  zodiacal  emis¬ 
sion  model  (ZEM)  is  also  needed.  Such  calibrations  and  a  ZEM  model  are  provided 
by  GEISHA. 


1  IRAS  DATA  ACCES.^ 

The  basic  saterl.  .  >r  the  GEISHA  project  is  the  raw  IRAS  data,  containing  all 
telemetry  data  fi'v  .  e  survey  instrument  detectors,  and  the  pointing  informa¬ 
tion.  The  dat'  were  iginally  taken  in  almost  random  scans  over  the  sky,  and 
presented  in  ordv  of  observation,  distributed  over  1000  tapes.  To  optimise  the 
GEISHA  system  for  detailed  investigation  of  (not  too  large)  sky  regions  both  the 
detector  data  for  the  scans  and  all  auxiliary  data  necessary  for  pointing  and 


86 


P  H  Wesselius  Ti  R  Boniekoe.  A  R  W  de  Jongt  and  DIM  Kesier 


flux  caliDralion  were  reorganises  into  a  system  of  363  plates.  A  plate  is  here 
defines  as  a  conerent  region  of  the  sky,  large  enough  to  fill  one  magnetic  tape 
in  FITS  format,  wnich  then  forms  a  totally  self-contained  dataset. 

For  routine  inspe''rion  of  the  data  a  set  of  programs  is  available,  applying  po¬ 
sition  ana/or  flux  calibration  to  the  raw  detector  outputs,  with  some  variation 
in  algorithm  uses.  The  end  products  of  these  programs  are 

-  a  file  ("sample  file")  containing  the  sampled  output  of  all  IRAS  detectors 
covering  the  region  of  interest  during  survey  scans,  each  sample  labeled 
with  time,  position  and  orientation  of  the  detector,  and  errors  in  these 
quantities ; 

-  an  image  of  the  region,  formed  by  the  surface  brightness  on  an  almost  arbi¬ 
trary  pixel  grid,  as  reconstructed  from  these  samples. 

All  intermediate  and  end  products  of  these  programs  are  in  FITS  format  to  ful¬ 
fill  Our  export  requirements.  The  subroutines  used  in  these  standard  programs 
are  available  as  a  software  library  for  special  research  projects. 

Currently,  we  have  the  detector  data  distributed  over  the  plate  tapes,  but  not 
yet  the  auxiliary  data.  Auxiliary  data  and  software  are  available  for  a  flux 
calibration  method  based  upon  tne  zodiacal  light,  positional  calibration  accors- 
ing  to  IPAC  methods,  and  imaging  by  co-adCing. 


II  EXPORTASILITl 


Exportability  was  formulated  as  a  secondary  design  goal  of  the  project,  both  for 
tne  data  sets  and  for  the  access  software.  This  would  give  the  system  better 
survival  chances  after  the  initial  effort,  enabling  us  to  migrate  to  the  best 
environment  for  f ur.ding.'computing/access  facilities. 

Tne  exportablity  was  achieved  by; 

-  use  of  Fortran-7V  as  the  commonest  language  in  astronomical  computing, 

-  isolation  of  computer- system  dependencies  into  a  small  set  of  subroutines, 

-  use  of  magnetic  tape  in  FITS-istandard  format  for  all  data  bases, 

-  use  of  two  widely  different  computer  systems  in  development,  a  DEC  PDP- 
11 /111  running  UNIX  V7  and  a  CDC  Cyber  170  running  NOS, 

-  a  (present)  transfer  of  the  GEISHA  user  system  to  a  VAX  8600  running  VMS. 
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III  CALIBRATION 

111.1  Introduce lor. 

The  GEISHA  system  will  supply  the  astronomer  with  a  choice  of  calibration  rou-' 
tines  and  zodiacal  emission  models  (ZEMs)  to  subtract  if  required. 

All  calibrations  have  a  common  structure:  it  is  assumed  that  the  raw  detector 
data  consist  of  a  baseline  {  the  detector  bias  )  plus  the  'true'  infrared  flux 
multiplied  by  a  gain.  Both  the  baseline  and  the  gain  might  vary  slowly  due  to 
photon  ana  particle  Induced  memory  effects,  instrumental  degradation  etc. 

The  main  differences  between  IPAC  (IRAS  Processing  and  Analysis  Centre,  U.S.A. ) 
and  GEISHA  are  in  philosophy.  IPAC  uses  the  instrumental  response  of  the  inter¬ 
nal  calibrator  as  tneir  standard  candle  to  determine  the  gain  and  the  daily  ob¬ 
servation  of  a  piece  of  SKy  near  the  ecliptic  north  pole  (TFPR)  as  a  monitor  of 
tne  baseline,  GEISHA  has  an  external  standard  candle:  the  zodiacal  emission. 
It  tries  to  autocal Ibrate  the  IRAS  database. 

Tne  first  calibration  schemes,  both  at  IPAC  and  GEISHA,  were  quite  simple:  the 
cnanges  in  gain  and/or  baseline  were  slow  and  global.  The  later  schemes  involve 
some  model  on  the  behaviour  of  the  detectors.  Different  instrumental  models 
give  rise  to  differences  in  calibration  routines. 

There  exists  a  whole  range  of  possible  ZEMs,  from  a  simple  lower-envelope  cubic 
spline  fit  to  a  full-fledged  physical  model  with  densities,  emissi vities,  and 
temperatures  as  a  function  of  position  in  the  zodiacal  dust  cloud.  Here  again 
the  choice  is  to  the  astronomer. 

111. 2  GEISHA  calibration  nr  i 

The  model  for  the  first  GEISHA  calibration  consists  for  each  scan  of  a  linearly 
Changing  baseline  plus  a  constant  gain.  The  zodiacal  emission  is  the  standard 
candle:  data  -  baseline  *  gam  *  ZEM.  A  gain  and  a  baseline  is  found  for  each 
detector  in  every  scan.  The  ZEM  itself  is  derived  as  the  mean  lower  envelope  to 
all  scans. 

As  the  intensity  of  the  zodiacal  emission  changes  also  with  the  solar  aspect  an¬ 
gle,  the  ZEM  can  not  simply  be  used  as  a  standard  candle.  In  this  calibration 
procedure  our  assumption  was  that  the  ZEM  could  be  separated  into  Che  product  of 
two  functionals,  one  in  ecliptic  latitude,  ZE90(B),  and  one  in  solar  aspect  an¬ 
gle,  f(e).  The  latter  is  temporarily  Ba,-ried  to  the  gain.  This  product  will  be 
derived  from  a  lower  envelope  fit  of  the  scan  to  the  ZE90.  Afterward  divorce 
takes  place;  the  e-dependent  part  moves  Into  the  ZEM  »  2E90(B)  *  f(e)  and  the 
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remainder  is  tne  gain. 

Tnis  calibration  is  operational  since  begin  1987. 

III. 3  IPAC  calibration  nr  1 

In  tne  first  IPAC  calibration  every  scan  has  a  constant  baseline.  Monitoring 
tne  TFPR  (total  flux  prime  reference  field)  twice  a  day  leads  for  every  SOP  (  - 
Satellite  Operating  Plan  -  1/2  day  )  to  a  new  set  of  baselines,  one  for  eacn 
detector.  In  flight  it  was  discovered  that  the  baselines  at  60  and  100  yra  were 
affected  by  tne  bias  boost.  A  slow  decay  was  included  with  timescales  of  50 
minutes  typically.  The  gam  is  calculated  from  the  response  on  the  internal 
calibrator  source  at  the  begin  of  tne  scan  and  at  the  end.  The  gain  was  allowec 
to  change  linearly  during  the  scan. 

All  presently  available  standard  IPAC  products  are  based  on  this  calibration. 
Note:  this  calibration  will  NOT  be  implemented  in  the  GEISHA  system. 

III.li  GEISHA  calibration  nr  2 

As  a  result  of  the  first  GEISHA  calibration  a  ZEM  was  found  that  was  a  product 
of  two  functions,  one  in  ecliptic  latitude  and  one  in  solar  aspect  angle.  The 
function  in  ecliptic  latitude  could  be  very  precisely  approximated  by  ar.  el-, 
lipse.  This  led  to  tne  conjecture  that  the  zodiacal  emission  as  a  whole  could  be 
approximated  by  the  emission  of  a  tilted  oblate  ellipsoid  of  uniform  density  and 
temperature.  At  the  IRAS  conference  in  London  (1987)  our  geometric  model  was 
corroborated  by  a  paper  presented  by  W.  Reach  and  C.  Heiles  who  had  made  a  phy¬ 
sical  ZEM  consisting  of  ellipsoidal  equidensl ties  and  had  the  model  integrated 
up  to  a  'last  equidensity  contour'. 

An  ellipsoidal  model  can  never  fit  the  so-called  sidebands  of  tne  zodiacal  emis¬ 
sion  (Hauser  et  al . ,  198*1).  The  sidebands  nave  been  identified  with  3  asteroid 
families  (Permott  et  al . ,  19811).  From  considerations  of  celestial  mechanics  it 
seems  reasonable  to  model  them  with  a  number  of  inclined  tori  fanned  over  ail 
nodal  angles.  A  side  band  model  (SBM)  Incorporating  these  ideas  is  still  under 
study  by  us. 

Apart  from  the  ZEM  there  is  always  a  contribution  of  galactic  emission.  The 
galactic  emission  is  so  granulated  that  only  a  map  of  the  whole  sky  can  function 
as  a  model.  A  fairly  coarse  map  is  sufficient.  We  took  an  equal  area  projection 
with  pixels  of  1  by  1  degree.  Each  pixel  is  the  mean  of  all  overlapping  detector 
values. 
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A3  a  standard  candle  «e  take  the  sum:  ZEM  ■*  SBM  *  GEM.  We  allow  the  baseline  to 
change  linearly  over  a  scan.  We  assume  that  the  slow  exponential  decay  as  found 
by  IPAC  will  De  sufficiently  approximated  by  a  straight  line.  The  non-linear 
behaviour  of  the  gam  will  be  corrected  according  to  a  model  developed  in  the 
report  "IR  Photon  Detectors"  May  7,  1986  by  Do  Kester. 

The  complete  model,  data  -  baseline  ♦  gain  •  (  ZEM  ♦  SBM  ♦  GEM  ),  contains  as 
free  parameters:  5  for  ZEM,  6  for  SBM,  40000  (pixels)  for  GEM,  4  for  each  scan 
(2  baseline  and  2  starting  conditions),  and  4  for  each  detector  (2  decaytimes 
and  2  increments).  In  total  about  64000  parameters  have  to  be  estimated,  still 
a  tiny  fraction  of  the  complete  database  of  13  Gbyte. 

This  calibration  and  the  pertaining  ZEM  will  be  available  by  middle  1988. 

111. 5  IPAC  calibration  nr  2 

In  IPAC's  second  model  the  baseline  is  a  slowly  varying  function,  now  for  the 
two  snorter  wavelength  bands  too.  Its  value  at  the  time  of  the  so-called  bias 
boost  decays  on  a  typical  timescale  of  1000  seconds.  For  the  two  longer 
wavelengths  (oO  and  100  um)  a  pnoton  induced  enhancement  is  incorporated  in  tne 
gam.  When  the  detector  response  exceeds  a  preset  threshold  the  gain'  increases 
proportional  to  the  detector  response.  That  increase  decays  on  timescales  of 
500  and  1500  seconds  for  50  and  100  ym  reap.  The  maximum  correction  is  10  to  20 
t. 

This  calibration  is  under  development  and  will  also  be  made  available  to  GEISHA 
users  by  middle  1988. 

A  physical  ZEM  is  under  development  by  J.  Good. 
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IV  SAMPLE  FILES 

In  order  to  make  the  raw  IRAS  data  in  a  computer-independent  way  available  to 
the  astronomical  community  a  new  product  has  been  made:  sample  flies.  For  a 
given  region  of  the  sky,  a  sample  file  contains  all  information  present  for  one 
of  the  IRAS  frequency  bands.  Each  detector  readout  is  present,  leading  to  a  long 
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taDie  Of  coorQinaces  and  brigncnesaes.  Tnis  table  naa  to  be  proceased  furtner  to 
maKe  maps  of  tne  region  of  interest.  Presently,  coaddea  maps  can  be  made  naving 
a  resolution  equal  to  tne  sizes  of  the  detectors;  about  5  arcalnute  in  tne 
cross-scan  direction,  and  0.75  to  3  arcminute  in  tne  in-scan  direction,  depend¬ 
ing  on  tne  frequency  band.  Tne  overall  accuracy  of  tne  coordinates  is  better 
tnan  3  arcseconds.  The  position  of  known  point  sources  is  reconstructed  only  to 
3"  in-scan  and  20"  cross-scan,  due  to  tne  asymmetrical  point  spread  function  of 
tne  detectors. 

Tne  brightness  can  be  calibrated  against  tne  zodiacal  light  emission  (GEISHA 
calibration).  Samplefiles  can  be  used  by  an  astronomer  who  wants  to  make  nis  own 
(nign-resQlution)  image  of  tne  sky. 


V  IMAGING 

V.i  Linear  Resolution  Enhancement 

In  oncer  to  give  an  Icea  of  the  amount  of  data  IRAS  gave  us:  1  square  degree  of 
tne  sky  near  tne  ecliptic  equator  has  been  sampled  i^OOOO  times  in  12  and  25  uit 
(oC  ym;  2COOO;  100  ym;  10000).  Towards  tne  ecliptic  poles  tnese  numbers  in¬ 
crease  by  an  order  of  magnitude,  due  to  tne  convergence  of  tne  scans. 

A  resampling  of  tne  AOOOO  detections  evenly  over  the  area,  would  give  a  linear 
sampling  distance  of  30  arcseconds.  A  linear  least  squares  reconstruction  tech¬ 
nique  has  been  devised  which  actually  yields  30  arcsecond  pixel  Images.  However, 
the  resolution  cannot  be  better  than  the  telescope  quality  permits,  i.e.  25", 
25",  60",  100"  at  12,  25,  60 ,  and  100  ytn.  The  resolution  will  oe  worse  because 
only  some  nigh  spatial  frequencies  are  present,  out  a  resolution  of  3  arcainutes 
Should  easily  be  attainable. 

The  method  works  as  follows.  If  tne  area  of  interest  is  divided  into  small  pix¬ 
els,  each  detector  readout  covers  a  number  of  them,  say  10.  The  measured  bright¬ 
ness  Is  tnen  tne  sum  of  the  unKnown  brightnesses  of  tnese  10  pixels  convolved 
with  the  point  spread  function.  This  coverage  can  be  described  as  one  equation 
in  many  unknowns  (-  all  pixels),  for  which  the  right-hand  side  is  the  one  meas¬ 
ured  brightness.  But  the  vast  majority  of  these  unknowns  have  a  coefficient 
equal  to  zero,  because  most  pixels  are  not  seen  by  the  detector  under  considera¬ 
tion.  The  next  readout  of  the  same  detector  now  covers  say  7  pixels  of  the  pre¬ 
vious  10,  and  3  new  ones  as  the  satellite  scans  the  sky.  This  new  equation  is 
partly  coupled  to  the  previous  one  through  the  7  common  unknown  brightnesses 
(which  are  assumed  to  be  constant). 
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The  coverages  of  the  area  should  Interconnect  all  pixels,  and  the  set  of  sparse 
equations  thus  formed  can  be  solved  for  the  unknown  pixel  values,  in  a  least 
squares  sense.  A  nessessary  condition  is  that  there  are  more  equations  (- 
readouts)  tnan  unknowns  (-  pixels).  Therefore,  in  principle  the  number  of  pixels 
per  unit  area,  and  thus  tne  resolution,  will  be  higher  near  the  ecliptic  poles 
than  at  tne  equator. 

V.2  Max..mum  Entropy 

A  further  enhancement  of  the  resolution  is  expected  from  Maximum  Entropy  Method 
(MEM)  image  restoration  techniques.  This  method  is  still  under  study,  but  a 
further  increase  of  the  resolution  to  30"  seems  possible  (Skilling,  private  com¬ 
munication).  The  resolution,  however,  need  not  be  constant  over  the  map:  areas 
with  more  coverages  will  be  reconstructed  with  higher  resolution  than  sparsely 
covered  ones.  A  second  map  of  the  final  resolution  then  has  to  be  provided. 

The  advantage  of  MEM  is  that  simultaneously  a  number  of  calibration  problems  can 
be  solved  for.  Presently,  when  a  strong  source  has  been  passed,  the  detector 
behaves  as  having  got  a  new  zero-point  and  a  new  sensitivity.  This  effect  is  not 
yet  corrected  for.  Also  slow  drifts  of  zero-.points  and  sensitivity,  apart  from 
the  ones  already  accounted  for  in  the  calibration,  can  cause  stripes  over  the 
map.  The  MEM  permits  parametrisation  of  these  effects,  and  solves  for  the  *best' 
image  and  the  'best'  determination  of  these  parameters  at  the  same  time.  The 
disadvantage,  however,  is  the  Increase  in  computational  cost  by  two  orders  of 
magnitude  for  the  same  area  of  sky. 


VI  CAPABILITIES  OF  GEISHA 

The  following  tasks  can  at  present  be  performed  by  the  user/astronoraer . 

-  Coordinate  transformations  between  any  of  the  following  systems:  Equatorial 
1950,  Equatorial  2000,  Galactic,  Ecliptic  1983.5,  and  Sunreferenced . 

-  Find  the  plate(s)  where  the  data  of  a  given  region  is  stored.  The  output 
is  one  or  more  platenumber (s) .  Snip  and/or  glue  a  new  plate  with  the 
desired  region  In  its  center.' 

-  Calibrate  the  scans  with  GEISHA  calibration  nr  1. 

~  Subtract  the  zodiacal  emission  node! . 

-  Plot  a  scan  as  a  function  of  time  or  of 

-  Perform  a  one-dimensional  coadd. 

-  Plot  the  run  of  all  scans  over  a  map. 

-  Coadd  all  scans  with  a  circular  beam  into  a  map. 

-  Make  a  sample  file,  i.e.  add  position  information. 
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-  Coadd  trie  sample  file  ultn  detector  sized/oriented  beam  Into  a  map. 

-  Make  a  contour  plot  of  a  map. 

wnenever  a  map  or  a  projection  is  mentioned  tne  astronomer  can  cnoose  Stereo- 
grapnic,  Cnoraonic,  Azimuthal  equidistant.  Azimuthal  equal  area,  Orthograpni c , 
Cylindrical  equal  area,  Mercator,  Sinusiodal  equal  area.  Hammer- Ai tof f  equal 
area  or  Cylindrical  equidistant. 

The  GEISHA  system  can  be  used  for  a  nominal  fee.  Call  31 "50-63^07^  and  asK  for 
Wesselius,  Kester  or  Arendz  to  arrange  a  visit. 
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1  Introduction 

The  Groningen  E.xpryrlabie  Infrared  Survey  High -resolution  Analysis  (GEISHA)  sys¬ 
tem  consists  of  a  collection  of  software  tools  which  allow  an  astronomer  to  create  in¬ 
frared  sky  maps  from  raw  IRAS  data.  The  Infra  Red  Astronomical  Satellite  (IRAS) 
observed  about  of  the  sky  by  a  semi-overlapping  scan  technique.  There  are  62 
detectors  in  the  focal  plane  in  4  wavelength  bands  at  12,  25.  60  and  100  pm.  (See 
IRAS  Explanatory  Supplement,  1985,  IR.AS  ES  hereafter) 

The  database  of  IR.AS  has  been  left  unaltered  apart  from  a  re-sort.  The  under¬ 
lying  idea  is  that  in  the  process  of  interpreting  a  certain  astronomical  object  new 
questions  may  arise,  which  require  a  different  processing  of  the  data.  Starting  from 
the  raw  data  has  already  proved  to  be  a  fruitful  strategy-. 

One  of  the  central  items  in  the  project  is  the  calibration.  The  procedure  for  the 
IR.AS  raw  data  calibration  consists  of  flatfielding.  For  each  detector  scan  a  set  of 
calibration  constants  is  extracted  by  fitting  the  lower  envelopes  of  the  detectors  of 
all  scans  to  each  other. 

2  IRAS  Data  Access 

The  basic  material  for  the  GEISHA  project  is  the  raw  IRAS  data,  containing  all 
preprocessed  telemetry  data  from  the  survey  instrument  detectors,  and  the  pointing 
information.  The  database  of  IRAS  contains  about  6000  scans  running  from  one 
ecliptic  pole  to  the  other  at  a  constant  angle  to  the  sun,  the  solar  aspect  angle. 
They  were  taken  at  an  angular  speed  of  3.75'/sec  and  at  a  rate  of  16,  16,  8  and  4 
times  per  second,  for  the  band  at  12,  25,  60  and  100  pm,  respectively.  The  scan 
strategy  was  organized  in  such  a  way  that  all  sources  were  seen  at  least  8  times  by  a 
detector,  yielding  a  data  volume  of  more  than  1000  tapes.  To  optimise  the  GEISHA 
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system  for  detailed  investigation  of  (not  too  large)  sky  regions  both  the  detector  data 
for  the  scans  and  all  auxiliary  data  necessary  for  pointing  and  flux  calibration  were 
reorganised  into  a  system  of  363  plates,  which  cover  the  whole  sky.  A  plate  is  here 
defined  as  a  coherent  region  of  the  sky,  large  enough  to  fill  one  magnetic  tape  in  FITS 
format  (Wells  et  al.,  1981),  which  then  forms  a  totally  self-contained  dataset. 

Currently,  we  have  the  detector  data  distributed  over  the  plate  tapes,  but  not  yet 
the  auxiliary  data.  Auxiliary  data  and  software  are  available  for  a  flux  calibration 
method  based  upon  the  zodiacal  light,  positional  calibration  and  imaging  by  co¬ 
addition.  The  positional  calibration  was  developed  at  the  IRAS  Processing  and 
Analysis  Center  (IPAC)  of  Caltech,  Pasadena,  USA.  (See  IRAS-ES) 

3  Calibration 

An  unavoidable  part  of  the  brightness  which  falls  on  the  detectors  is  the  zodiacal 
emission.  It  originates  from  dust  inside  the  solar  system  and  provides  a  broad  back¬ 
ground  (or  better  foreground)  feature  with  a  FWHM  of  »  50°.  It  is  dominant  at  12, 
25,  and  60  pim,  and  somewhat  less  at  100  fim.  Although  it  is  a  nuisance  for  most 
astronomers,  it  can  be  used  as  a  standard  candle  for  calibration  purposes.  In  the 
same  process  a  zodiacal  emission  model  (ZEM)  can  be  established. 

The  GEISHA  system  will  supply  the  astronomer  with  a  choice  of  calibration  rou¬ 
tines  and  zodiacal  emission  models  (ZEMs)  to  subtract  if  desired.  All  calibrations 
have  a  common  structure:  it  is  assumed  that  the  raw  detector  data  consist  of  a  de¬ 
tector  bias,  plus  the  ‘true’  infrared  brightness  multiplied  by  a  gain  factor  g.  Both 
the  bias  and  the  gain  might  vary  slowly  due  to  photon  and  particle  induced  memory 
effects,  instrumental  degradation  etc.  A  whole  range  of  possible  ZEMs  can  be  devel¬ 
oped,  from  a  simple  lower-envelope  cubic  spline  fit  to  a  full-fledged  physical  model 
with  densities,  emissivities,  and  temperatures  as  a  function  of  position  in  the  zodiacal 
dust  cloud.  Here  again  the  choice  is  to  the  astronomer. 

The  total  IRAS  database  is  clearly  too  large  to  handle  in  one  program.  The  first 
step  in  the  reduction  of  the  database  was  a  cubic  spline  fit,  which  brought  the  sample 
rate  down  to  one  datapoint  per  second  (or  3.75  arcrain).  The  splines  are  part  of  the 
GEISHA  system;  they  can  be  used  for  the  study  of  large  scale  structures.  In  toted  it 
is  still  1  Gbyte. 

The  second  step  which  was  only  taken  for  calibration  purposes,  involved  a  lower 
envelope  fit  over  one  degree  to  the  splines.  The  lower  envelope  fit  entailed  that  all 
single  point  sources  dropped  out;  we  are  mainly  interested  in  the  large  scale  structure 
of  the  zodiacal  emission.  This  last  database  contains  70  Mbyte  or  18  Mbyte  in  each 
of  the  4  bands. 

4  Linear  Model 

In  the  first  GEISHA  calibration  a  linearly  changing  bias  plus  a  constant  gain  was 
assumed  for  each  detector  scan. 
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O  =  a  +  6<+ (?(ZEM  +  SKY),  (1; 

where  O  is  the  detector  output,  a  +  bt  is  the  linearly  changing  bias  (the  baseline), 
g  is  the  gain  factor,  ZEM  is  the  zodiacal  emission  and  SKY  is  the  astronomical  sky. 
As  the  astronomical  sky  is  quite  empty  in  the  infrared,  the  SKY  term  will  be  ignored 
in  the  linear  model.  The  zodiacal  emission  is  a  function  of  the  position  of  the  earth, 
and  the  pointing  direction  of  the  satellite,  given  by  the  solar  aspect  angle  0,  and  the 
angle  along  the  scan  direction,  The  angle  6  is  constant  during  a  scan.  In  a  first 
approximation  the  annual  motion  of  the  earth  through  the  somewhat  tilted  dustcloud 
(Hauser  et  al.,  1984)  can  be  corrected  with  a  small  shift  in  rl>. 

For  all  positions  on  a  scan  the  brightness  of  the  zodiacal  emission  increases 
with  decreasing  solar  aspect  angle;  closer  to  the  sun  the  zodiacal  emission  becomes 
brighter.  It  is  assumed  that  the  ZEM  can  be  written  as  the  product  of  two  functions: 

ZEM  = /(fl)ZEM9o(v'>),  (2) 

where  ZEMgoCV'')  is  the  zodiacal  emission  model  at  a  solar  aspect  angle  ff  =  90°,  and 
f(9)  is  a  scaling  function.  Both  functions  are  as  yet  unknown  and  will  be  estimated 
from  the  data  by  an  iterative  procedure. 

Knowing  the  n-th  iterate  of  the  function  ZEM9o(v)n  we  observe  that  for  each 
individual  detector  scan  the  output  can  be  written  as  a  linear  equation  in  the  3  pa¬ 
rameters  which  have  to  be  estimated:  a„+i,  6„+i  and  the  combination  of  {g  fi6))n+i- 

O  =  a„+i  -f  bn+\  t  +  (g  f(9))n+l  ZE.M9o(l'*)n.  (3) 


All  thus  calibrated  detector  scans  are  averaged  into  the  next  iterate  of  ZEM9o(t/')n+i. 


ZEM9o(t!')n-(-l 


^  ^n+ 1  ^n+ 1 ^ 


(P /(»)).+  ! 


(4) 


In  both  fitting  procedures  we  are  only  interested  in  the  common  lower  envelope. 
All  datapoints  which  stray  from  this  lower  envelope  we  consider  as  spurious.  They 
are  eliminated  with  an  object  function  which  is  more  robust  than  the  normail  least 
squared  method.  We  took  a  one-sided  biweight  function  of  Tukey  (Goodal,  1983). 

When  it  is  sufficiently  converged,  typically  after  about  10  iterations,  the  factor 
g  f{0)  is  separated  for  each  detector  into  a  part  f{9)  which  comprises  the  common 
functional  behaviour  with  6,  and  a  remaining  part.  Th’  former  is  multiplied  with 
ZEM90  and  represents  the  ZEM;  the  latter  is  the  g2un  for  each  detector  scan. 

This  calibration  is  operational  since  begin  1987. 


5  Detector  Behaviour 

IR  photon  detectors  of  the  kind  IRAS  was  supplied  with  are  known  to  show  particle 
and  photon  induced  memory  effects  (Kester,  1986).  A  full  physical  model  of  these 
effects  still  eludes  all  efforts  (Westervelt  and  Teitsworth,  1985).  The  model  here 
presented  only  concerns  photon  induced  memory  effects  on  the  gain  factor.  The 
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particle  induced  effects  due  to  the  passage  through  the  South  Atlantic  Anomaly 
seem  to  influence  the  bias  mainly  (See  IRAS-ES).  We  already  allow  for  a  linearly 
changing  bias  and  hope  that  that  will  suffice. 

The  non-linear  behaviour  of  the  gain  will  be  corrected  according  to  a  two  energy 
level  model  for  IR  photon  detectors.  The  model  can  be  described  as  two  coupled  first 
order  differential  equations. 

cr~  +  h  A  F  hm,n.  (5) 

at 


and 


T-^  +  g  =  B  F  +  gmxn  (6) 

The  equations  are  coupled  by 


Here.  g{t)  is  a  time  dependent  gain  which  couples  the  detector  output  O  to  the  total 
incident  radiation  F  according  to  equation  (1).  Further,  h{t)  is  a  function  which  is 
inversely  proportional  to  the  decay  time  r(t)  of  the  gain.  It  acts  as  a  saturation  on 
the  gain  at  large  input  fluxes.  The  decay  time  of  the  saturation  function,  hit),  itself 
is  (T.  The  constan's  a,  tq,  A,  B,  hmm  and  g^tn  are  parameters  to  be  estimated:  one 
of  them  can  be  chosen  freely.  We  take  =  1-  1  he  fact  that  h  is  always  greater 
than  or  equal  to  1  entails  that  g  changes  on  time  scales  equal  or  smaller  than  r^. 

To  estimate  the  parameters  of  the  gain  model  the  two  differential  equations  must 
be  rewritten  as  difference  equations.  (Ljung,  1987). 


hk  =  (I  -  p){AFk-i  +  \)  A  ph^^i,  (8) 

9k  =  [I  -  g){BFk-i  A  grn.rx)  A  q  9k-}.  (9) 

w  here 

p  =  exp{-T/c),  (10) 

q  =  exp{-Thk-}/To).  (11) 


T  is  the  sample  time  interval.  From  these  equations  an  estimator  of  the  momentary 
gain,  Qk,  is  found  in  terms  of  the  previous  /i*_i  and  gk-j,  and  in  terms  of  the  previous 
input  fluxes  Fk-i  and  Fk-2.  Via  equation  (1)  an  estimator  of  O*  is  found  which  can 
be  compared  with  the  actual  detector  scan  values.  The  difference  is  to  be  minimized 
under  a  robust  objective  function. 
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For  most  of  the  scans  the  linear  scheme  proved  to  be  reasonably  successful.  But 
there  are  broken  scans,  running  only  part  of  a  semicircle,  for  which  the  baseline  and 
the  gain  can  not  be  separated  any  more.  This  sometimes  results  in  negative  gains 
or  other  kinds  of  bad  fits,  particularly  when  the  short  scan  covers  some  remaining 
gtilactic  emission.  Moreover  for  scans  taken  at  extreme  solar  aspect  angles  [6  less 
than  75°  or  more  than  105°),  the  assumption  of  separability  of  ZEM,  as  stated  in 
equation  (2),  does  not  hold. 

Still  a  zodiacaJ  emission  model  was  found  of  which  the  ZEM9o(ti’)  pa-rt  could 
be  very  precisely  approximated  by  an  ellip.s*'.  This  led  to  the  conjecture  that  the 
zodiacal  emission  as  a  whole  could  be  approximated  by  the  emission  of  a  tilted  oblate 
ellipsoidal  dust  cloud  of  uniform  density  and  temperature.  The  ZEM  is  a  function  of 
solar  aspect  angle  9,  the  longitude  of  the  sun-referenced  coordinate  system  t/’,  and  the 
position  of  the  earth  represented  by  the  solar  elongation  This  idea  is  supported 
by  Reach  and  Heiles  (19S7). 

An  ellipsoidal  model  can  never  fit  the  so  called  sidebands  of  the  zodiacal  emission 
(Hauser  et  al.,  19.81).  The  sidebands  have  been  identified  with  3  asteroid  families 
(Dermott  et  al.,  198-1)  From  considerations  of  celestial  mechanics  (Hirayama.  1918 
and  Brouwer,  1951)  it  seems  reasonable  to  model  them  with  a  number  of  inclined 
ton,  fanned  over  all  nodal  angles.  A  side  band  model  (SBM)  incorporating  these 
ideas  is  still  under  study. 

Apart  from  the  zodiacal  emission  there  is  always  a  contribution  of  astronomical 
background  (SK"^');  most  of  it  is  galactic  emission.  The  aistronomicai  background  is 
so  granulated  that  only  a  map  of  the  whole  sky  can  function  as  a  model.  A  fairly 
coarse  map  in  an  equad  area  projection  with  pixels  of  1°  by  1°  will  suffice.  Each  pixel 
is  the  mean  of  all  overlapping  detector  values.  As  a  standard  candle  the  sum  of  ZEM, 
SBM  and  SKY  can  be  taken,  while  g{t)  follows  the  detector  model  expounded  in  the 
section  5. 

The  complete  mod'll, 

0  =  a  +  6/  +  (7(0(ZEM(v\tf,A^J  +  SBM(t/-,^,  Ag)  +  SKY(A, /?)),  (12) 

contains  as  free  parameters:  5  for  ZEM,  6  for  SBM,  40000  pixels  for  SKY,  5  for  each 
scan  (  2  baseline  and  3  starting  conditions  ),  and  5  for  each  detector  (  2  decaytimes, 
2  increments  and  a  minimum  gain  ).  In  total  about  70000  parameters  have  to  be 
estimated,  still  a  tiny  fraction  of  the  complete  database  of  13  Gbyte. 

The  non-bnear  model  asks  for  a  more  compbeated  iterative  scheme.  Each  iteration 
requires  three  passes  over  the  data:  one  to  determine  SKY,  one  for  ZEM  +  SBM  and 
one  for  the  baseline  and  the  gain.  In  each  pass  the  items  which  are  not  fitted,  attain 
the  values  derived  in  the  previous  iteration.  SKY  is  found  by  coadding  the  detector 
scans,  properly  scaled  with  baseline  and  gain  and  stripped  of  the  contribution  of 
ZEM  +  SBM,  into  an  equal  area  map  of  the  whole  sky.  ZEM  and  SBM  are  non- 
bnear  functions  of  solar  elongation,  sun-referenced  longitude  and  latitude.  An  elegant 
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way  to  fit  non-linear  funoticns  is  the  Leve'^berg-Marquardt  me'hod  (e.g.  Press  et  aJ. 
1986). 

The  derivation  of  the  calibration  paiaineters  falls  into  two  parts.  For  each  detector 
the  values  of  the  gairs  model  have  to  be  estimated  in  the  way  that  was  outlined  in 
the  sect'ion  5.  And  for  each  detector  scan  the  baseline  has  to  be  estimated  plus 
the  fictitious  values  b  ,  go  and  fo.  which  represent  the  state  of  the  system  at  the 
beginning  of  the  detector  scan. 

Most  ot  the  steps  outlined  above  need  iterations  in  themsehes.  E.g.  c*.  nor-linear 
regression  model  cm  only  be  solved  i.i  ar.  iterative  way  Nontheless  it  is  hoped  that 
one  grand  iteration  loop,  in  winch  each  set  of  parameters  is  estimated  only  once 
during  an  iteration  loop  will  eventually  converge.  The  choice  of  starting  values  fo* 
the  parameters  may  be  quite  im5:>ori<i.nt. 

Thus  calibr-i’ion  and  the  jierlaining  ZE.M  will  be  available  by  middle  1988. 

7  Further  Tools 

For  routine  inspection  of  the  data  a  set  of  programs  is  acailable.  applying  position 
a::i;'. /or  fiu.x  calibration  to  the  raw  detector  outputs,  with  some  variation  in  algorithm 
u.-ed.  1  he  <  r.d  products  of  these  programs  are 

•  a  fib:  i  "sample  file'  i  containing  the  sampled  output  of  all  IRAS  detector.'  cov¬ 
ering  the  region  of  uitert,,st  during  survey  scans,  each  sample  labeled  wiili  lime, 
po'Uioi.  and  orientation  of  the  detector,  and  errors  in  these  qua.itities; 

•  moderate  resolution  images,  formed  by  averaging  of  the  overlapping  lamples 
on  an  almost  arbitrary  pixel  grid; 

•  high  resolution  images  obtaine.l  by  a  linear  least  squares  reconstruction  tech¬ 
nique  exploiting  the  overlap  between  the  sample.',  and  the  redundancy  in  the 
co'.orage  of  the  'ky. 

Ail  ii; lerinediate  and  end  products  of  these  programs  are  in  FITS  formal  to  ful¬ 
fil'  our  export  requirements.  The  subroutines  used  in  these  standard  programs  are 
avaiiah'''-  a:-  a  .-•'ftware  library  for  special  research  projects. 
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Image  Reconstruction  from  IRAS  scans 
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Captions 

•  Caption  Figurf  1:  IRAS  focal  plane.  Of  the  62  infrared  detectors  the  3  filled-in 
detectors  were  inoperative. 

•  Caption  Figure  2:  Point  spread  function  of  detector  15  (60  pm).  The  contour  levels 
are  at  99,  98,  95,  90  -  10  (10),  5,  2,  and  1  percent  of  the  maximum.  The  nominal  size 
is  represented  by  the  rectangle,  which  corresponds  to  the  rectangle  in  Figure  1. 

•  Caption  Figure  3:  Detector  outputs  of  one  scan  covering  M51  (60  ^m).  The  scan 
direction  is  from  right  to  left,  and  the  scan  has  a  length  of  one  degree.  The  detector 
outputs  arc  displaced  vertically  corresponding  to  their  cross-scan  position  in  the  focal 
plane.  The  figure  is  corrected  for  the  relative  positions  of  the  detectors  in  the  in-scan 
direction  and  can  be  regarded  as  a  ruled  surface  plot.  The  brightnesses  have  been 
scaled  relathe  to  the  maximum  (set  at  100)  in  this  piece  of  the  scan. 

•  Caption  Figure  4'  Positions  of  the  samples  in  the  area  of  M51  (60  pm).  The  circles 
represent  the  two  smallest  detectors  in  the  band,  viz.  detectors  11  and  31  (see  Fig.  1). 
The  scans  roughly  run  from  top-left  to  bottom-right.  The  area  is  covered  by  about 
700  samples,  and  the  nominal  size  of  one  detector  is  outbned.  A  grid  with  a  one 
arcminute  spacing,  suitable  for  image  reconstruction,  is  superimposed.  The  positions 
are  projected  on  a  plane  tangent  to  the  map  centre. 

•  Caption  Figure  5:  A  ‘co-added’  map  of  M51,  i.e.  using  the  (5x1  arc  minutes)  detector 
size  at  60  pm  as  the  beam.  The  size  of  the  map  is  15  by  15  arc  minutes. 

•  Caption  Figures  6:  Maps  of  M51  (60  pm)  using  the  linear  least-squares  reconstruction 
method  with  two  different  damping  factors. 

•  Caption  Figure  7:  A  CPC  image  of  M51  (50  pm).  The  resolution  of  this  map  is 
1 .5  X  1 .5  arc  minutes. 

•  Caption  Figure  8:  A  ‘FIER’  reconstruction  of  an  M  51  AO  observation  (IRAS,  60  pm). 
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LINEAR  LSQ  IMAGE  RECONSTRUCTION  M51  AT  60ptn 
"DAMPING"  =0.10 
FLUX  IN  MJy/SR 


IRAS  CPC  IMAGE  OF  M51  AT  50pm 
FLUX  IN  MJy/SR 
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A  maximum  entropy  "FIER"  restoration  of  the  60pm 
IRAS  AO  data  of  M51,  done  by  MRC/TUFTS. 
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